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By means of the recently proposed algorithm based on the tensor product states, the magnetization 
process of the spin-1/2 anti-ferromagnetic XXZ model on a square lattice is investigated. In 
the large spin-anisotropy limit, clear evidence of a first-order spin-flip transition is observed as an 
external magnetic field is increased. Our findings of the critical field and the discrete jumps in 
various local order parameters are in good agreement with the quantum Monte Carlo data in the 
literature. Our results imply that this algorithm can be an accurate and efficient numerical approach 
in studying first-order quantum phase transitions in two dimensions. 
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Numerical simulations are usually required in the the- 
oretical investigation on strongly correlated systems, be- 
cause analytical solutions are not available in most cases. 
Consequently, developing accurate and efficient numeri- 
cal tools becomes one of the central issues in the un- 
derstanding of quantum many-body systems. Recently, 
based on an efficient representation of two-dimensional 
system's wave function through a tensor network, a se- 
ries of new simulation algorithms has been achieved. In 
particular, the infinite projected entangled-pair states 
(iPEPS) algorithm^ has been proposed and applied to 
various interesting systems with success In this 

approach, the ground-state wave function is described 
by the so-called tensor product state (TPS)- ,§ or the 
projected entangled-pair state (PEPS)i2iiS Taking into 
account possible translational symmetry in the ground 
state, such a tensor network can be simply represented 
by copies of a small number of tensors even for systems 
on infinite lattices. After optimizing these tensors under 
specific prescriptions, a number of physical properties can 
be calculated from the optimized TPS/PEPS. 

By handling tensor-product wave functions in different 
manners, schemes distinct from iPEPS algorithm have 
also been put forward. n i 12 A virtue of these approaches 
is that they can be implemented with ease. In Ref. [Til , 
the optimized TPSs are determined via direct variational 
approach, where the variational energies of systems of 
very large sizes are efficiently evaluated by means of the 
tensor renormalization group (TRG) method. 13 ' 14 The 
expectation values of physical quantities are then cal- 
culated from the optimized TPS again under the TRG 
method. This algorithm has been tested for several two- 
dimensional (2D) quantum spin models^ and the re- 
sults agree well with previous findings. Alternatively in 
Ref. Qjl, the ground states of a TPS form are obtained by 
using the power method through iterative projections. 
This approach can be considered as a generalization of 
the ID infinite time-evolving block decimation (iTEBD) 
methodic to the two dimensional cases. After getting the 
ground states, the TRG method^ 3 - is employed to calcu- 



late the expectation values of physical observables. It is 
shown that accurate results for the Heisenberg model on a 
honeycomb lattice can be reached under this approach. 12 

Due to the simplicity and efficiency of the iTEBD and 
the TRG algorithms, the approach proposed in Ref. Q~l 
can become one of the promising numerical methods 
in studying quantum many-body systems once its gen- 
eral validity is established. Recently, it is shown that 
TPS/PEPS ansatz is suited to study the first-order phase 
transition. 3 However, because of the difference in opti- 
mizing ground states and in evaluating expectation val- 
ues, one may wonder if the combined iTEBD and TRG 
algorithm can determine the first-order phase transitions 
to the same accuracy as the iPEPS algorithm does. 

In order to provide further benchmark on the perfor- 
mance of the combined iTEBD and TRG algorithm, in 
this work we investigate the magnetization process of the 
spin-1/2 anti- ferromagnetic XXZ model on a square lat- 
tice. Here the large spin-anisotropy case is considered, 
where the existence of first-order spin-flip transitions in 
the magnetization process has been established by means 
of quantum Monte Carlo (QMC) simulationsJ&iZiiS We 
find that various local order parameters defined below 
change discontinuously at a critical field, which clearly in- 
dicates the appearance of a first-order transition. More- 
over, satisfactory results of the critical field and the dis- 
crete jumps in the local order parameters are be ob- 
tained as compared to the previous QMC findings.— Our 
present investigation suggests that this combined algo- 
rithm should also be an effective numerical method in 
studying first-order quantum phase transitions in two di- 
mensions. 

Before presenting our results, it is instructive to 
sketch the combined iTEBD and TRG algorithm em- 
ployed here. We know that the ground state can in 
principle be determined through the imaginary time 
evolution for a given initial state \^o): I^gs) = 
lim T ^ 00 exp(-Fr)|*o)/||exp(-ffr)|*o)||- If, just like 
the present case, the model Hamiltonian can be written 
as a sum of terms huj\ involving only pairs of nearest- 



neighboring sites i and j, the Suzuki- Trotter formula^ 
can be exploited to decompose the imaginary time evo- 
lution operation into a product of two-site evolution op- 
erators: = exp(— h(j j)St), where St <C 1. It is 
also known that any wave function can always be ap- 
proximated in a TPS form. A possible construction of 
TPS for systems on a square lattice is to attach a rank- 
five tensor [Fj]* rjtd to each site i and a diagonal singular 
value matrix (hence a vector) [A/jjv]/ to each bond of 
nearest- neighboring sites i and j. Here s is the physi- 
cal index with s = 1, 2 for the present spin- 1/2 case, and 
I, r, u, d{= 1 • • • D) denote the virtual bond indices in four 
directions. In general, better representation of a given 
wave function can be achieved by increasing the bond 
dimension D. Taking into account the possible transla- 
tional symmetry in the ground state under shifts by two 
lattice sites both in the x and y directions, the tensor 
network can be simply represented by copies of tensors 
within a 2 x 2 unit cell. That is, we are left with four inde- 
pendent Ti tensors and eight independent Xnj) matrices. 
The action of a two-site evolution operator Unj) on such 
a TPS can be absorbed by performing a singular value 
decomposition, and thus leads to an update of the Ti, 
Tj, and X^j\ tensors^ When eight nearest-neighboring 
bonds within the 2x2 unit cell are all updated, a com- 
plete iteration is achieved. After sufficient time of such 
updating iterations, the optimized ground state of the 
TPS form can be generated. 

Since evaluation of the expectation values for a TPS 
under the most straightforward method is exponentially 
difficult, for a complete numerical algorithm, an effi- 
cient way to do these calculations for large systems must 
be also constructed. Here the TRG approac h 11 ' 13 is 
employed. For any operator that can be decomposed 
into product of local operators, O = fl^Oi, evaluating 

(*gs|0|*qs) for the TPS ground state |^gs) is equiva- 
lent to compute the contraction of a corresponding tensor 
network of T tensors. Within such a tensor network, the 
rank-four tensor T l at site i is defined as 

[T l W = J2i S '\ 6 ^ A ^rud (1) 

ss' 

with 



[^Ifrud = y[\i,i-x)]l[^{i,i+x)]r[^(i,i+^)]u[\i,i-H)]d 

X [^illrud ; (2) 

where \s) represents the spin state at site i. i ± x and 
i ± y denote the nearest neighbors of site i in the x and y 
directions, respectively. I = (I, V) is the double bond in- 
dex, and f, u, d are similarly defined. The tensor network 
of T tensors then can be coarse-grained in an iterative 
fashion i 11 ' 13 Each complete renormalization group (RG) 
step reduces the size of the network by a factor of 2. The 
accuracy of such a RG process is controlled by a cutoff 
D cut on the double bond indices of the coarse-grained 
tensor. Therefore, to evaluate the contraction of the ten- 
sor network of size 2™ +1 x 2" +1 , we need only perform n 
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FIG. 1: (Color online) Values of mj, m£, and ml for the 
ground state at A = 1.5 as functions of external field h for 
systems of size 2 7 x 2 7 with D — 4 and D C ut = 16. 

RG steps. To sum up, the TPS provides an efficient way 
to approximate the 2D wave functions. The agreement 
between the actual wave function and the represented 
TPS wave function can be improved simply by increas- 
ing the bond dimension D. Besides, the TRG approach 
serves as an efficient tool to evaluate the expectation val- 
ues for a TPS ground state of very large systems, where 
the accuracy can be systematically improved by increas- 
ing the cutoff Z) cu t. In the present work we consider the 
bond dimension up to D = 5 and keep D cut > D 2 to 
ensure the accuracy of the TRG calculation. 

The general simulation procedure is described as fol- 
lows. For a given h and D, we take a set of random T 
and A tensors as our initial state \^o). While the ini- 
tial state may not have the spatial rotational symmetry, 
during the imaginary time evolution, the evolved state 
will converge towards a ground state which respects this 
expected symmetry. It hence provides a self-consistent 
stability check for the algorithm. To minimize the Trot- 
ter error, we usually start with St = 10 _1 and gradually 
decrease it to St — 10 -3 to ensure the convergence of the 
wave function. 

In the following, we present our numerical results for 
the spin-1/2 XX Z model with systems size N — 2 7 x 2 7 . 
In the presence of an external magnetic field h along z 
direction, the Hamiltonian of the XXZ model is given 
by22 

(id) i 

where Sf is the a{= x,y,z) component of the spin-1/2 
operator at site i, and (ij) runs over all the nearest- 
neighboring pairs of spins at sites i and j. J = 1 is the ex- 
change coupling, and A(> 0) is an anisotropic parameter. 
We focus our attention on the large spin-anisotropy case 
of A = 1.5, where accurate QMC calculations have been 
performed^ The expectation values of the ^-component 
staggered magnetization m z s = J2i(Sf)e lCi ' r '/4, the uni- 
form one m z u = ^2AS£) /£, and the x-component uniform 
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FIG. 2: (Color online) Energies e(h) per site for the adia- 
batically evolved states |*l/_t(/i)) and {^hQi)). The values for 
the ground state I^gsC*)) are denoted by open circles. Here 
A = 1.5 for systems of size 2 7 x 2 7 with D — 4 and D cut = 16. 
The inset shows the critical field h c for various D (D cu t = 16 
for D < 4 and D cu t = 25 for D = 5). Dotted line is guide to 
eyes. 



magnetization = ^2,ASf)/4: for the ground states 
I^gsC 1 )) are shown in Fig. [TJ where Q = (tt, 7f) and 
the sum on i runs over four sites within the 2x2 unit 
cell under consideration. Here we take the bond di- 
mension D = 4 and the TRG cutoff £> cut = 16. We 
note that results for D = 3, 5 are very similar to those 
for D = 4 and are thus not shown here. In the large 
spin-anisotropy limit with A > 1, it is known that a 
first-order spin-flip transition from a Neel-ordered phase 
to a spin-flopping phase will occur as the external field 
h increases from zeroj 16 > 17 i 18 Crossing the critical field 
h c , the z-component staggered magnetization m z sud- 
denly drops to zero and the uniform part m z u jumps to 
a nonzero value. When h > h c , m z u increases monotoni- 
cally and finally reaches its saturated value {m z u = 1/2) 
at h = h s [= 2(1 + A)], while the staggered one m z re- 
mains zero. A character of the spin-flopping states for 
h c < h < h s is the existence of finite x-componcnt mag- 
netization , whose value gives a measure of spin super- 
fluidity (see below). When spins become fully polarized 
in z direction as h approaches h s , will decrease to 
zero. As seen from Fig.[TJ our values of m z , to*, and mjj 
do show the expected results. 

Typically a first-order quantum phase transition comes 
from energy level crossing in the ground state, and the 
crossing point gives the critical value of the tuning pa- 
rameter^ In the present case, the relevant states should 
be the ground state in the Neel-ordered phase with zero 
m z u and that in the spin-flopping phase with finite toJ. 
Here we simulate the adiabatically evolved states ^^(h)) 
and \^n(h)) starting from the computed ground states 
in the Neel-ordered phase and in the spin- flopping phase, 
respectively. That is, \^z(h)) (\^ji(h))) are determined 
starting from the ground state I^gsCMiu)) for a given 
initial parameter h lnl < h c (/i; n i > h c ), and adiabatically 
increasing (decreasing) h in the Hamiltonian well beyond 
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FIG. 3: (Color online) Values of (a) m z u (b) \m z a \ 2 (c) \m x u \ 2 
for the adiabatically evolved states \^l{Ji)) and \^r{K)). The 
values for the ground state |^gs(/i)) are denoted by open 
circles. Here A = 1.5 for systems of size 2 7 x 2 7 with D = 4 
and Dcut = 16. 



crossing the critical field h c . True ground states | ^Gs(h)) 
are the ones with lower energies. The corresponding en- 
ergies e(h) per site are shown in Fig. [5] for D = 4 and 
£>cut = 16. We find that the energies of ^^(/i)) remain 
unchanged as h varies, while those of \^>u(h)} are lowered 
as h increases. This is expected since \^>L{h)) should be- 
have like Neel-ordered states with zero m z u and thus their 
energies do not depend on the external field. However, 
\^r(K)) should retain the spin-flopping character with 
nonzero to*, hence their energy expectation value (H) for 
the Hamiltonian in Eq. © should behave as a decreasing 
function of h. Due to level crossing in these two states, 
discontinuity in the first derivative of the ground state 
energy appears. This again indicates the presence of a 
first-order quantum phase transitions. From the crossing 
point in Fig. [2l we find that h c ~ 1.829, which is quite 
close to the value estimated by QMC (h c ~ 1.83). 1 — The 
dependence of h c on the bond dimension D is plotted in 
the insect of Fig. [5J While the findings of h c for D = 2 
and 3 are somewhat higher than the QMC result, satis- 
factory values can be obtained for larger D. 

To show further evidence of a first-order transition 
between the Neel-ordered phase and the spin-flopping 
phase, the results of m*, |m z | 2 , and |m^| 2 for the adi- 
abatically evolved states \\&L(h)) and \^n(h)), and the 
corresponding values for the ground state I^gsW) are 



displayed in Fig. [31 We find that the values of m z ul 
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and |mjj| 2 for the ground state are all discontinuous at 
h c . Moreover, all the results for the adiabatically evolved 
states and |\I/_r(/i)} show clearly the hysteresis 

behaviors. These facts strongly support the presence of 
a first-order transitions. The discrete jumps at h c for m z u 
and \m z | 2 are c ~ 0.125 and |mf iC | 2 ~ 0.202, respec- 
tively. Both of them agree with the QMC results reported 
in Ref. [l7|: m z uc ~ 0.11 and \m z sc \ 2 ~ 0.20^2. As men- 
tioned before, in the spin- flopping phase for h c < h < h s , 
there exists spin superfluidity which can be characterized 
by nonzero spin stiffness p s (or the superfluid density in 
the corresponding hard-core boson model). It is shown 
that p s changes discontinuously at the first-order spin- 
flip transition.— To the best of our knowledge there is 
no straightforward way to calculate the spin stiffness p s 
within the TPS framework. Instead, the square of the x- 
component magnetization are evaluated, which can 
be related to the density of Bose condensate in the corre- 
sponding hard-core boson models Thus the fact of non- 
vanishing \m^\ 2 also implies the existence of spin super- 
fluidity. As seen from Fig.[3Jc), similar to the behavior of 
p s observed in the QMC study, \m*\ 2 has also a discrete 



jump at h c , which is of magnitude J 2 = 0.118. 

In summary, the first-order spin-flip transition of the 
spin-1/2 XX Z model with large spin anisotropy can be 
detected under the combined iTEBD and TRG algorithm 
proposed in Ref. Il2j Good agreement with the accu- 
rate QMC calculations can be obtained by using merely 
moderate bond dimension D and the TRG cutoff D cut . 
This demonstrates that the current formalism will be a 
competitive numerical method to determine particularly 
first-order quantum phase transition in two dimensions, 
with the simplicity and efficiency as its advantage. We 
note, however, that further investigations are necessary 
to establish its general validity, and to explore its rela- 
tive performance as compared to other TPS/PEPS-based 
approaches. 

C.-Y. Lai and P. Chen thank the support from the 
National Science Council of Taiwan under Contract No. 
NSC 95-2112-M-007-029-MY3. M.-F.Yang acknowledges 
the support by the National Science Council of Taiwan 
under Grant No. NSC 96-21 12-M-029-004-MY3. This 
work is supported by NCTS of Taiwan. 



* Electronic address: pcchen@phys.nthu.edu.tw 

1 J. Jordan et ai, Phys. Rev. Lett. 101, 250602 (2008). 

2 H.-Q. Zhou, R. Orris, and G. Vidal, Phys. Rev. Lett. 100, 
080602 (2008). 

3 R. Orris, A. Doherty, and G. Vidal, Phys. Rev. Lett. 102, 
077203 (2009). 

4 B. Li, S.-H. Li, and H.-Q. Zhou, Phys. Rev. E 79, 
060101 (R) (2009). 

5 J. Jordan, R. Orris, and G. Vidal, Phys. Rev. B 79, 174515 
(2009). 

6 B. Bauer, G. Vidal, and M. Troyer, J. Stat. Mech.: Theory 
Exp. (2009) P09006. 

7 M. A. Martm-Delgado and G. Sierra in Density-Matrix 
Renormalization - A New Numerical Method in Physics, 
Lecture Notes in Physics Vol. 528, edited by I. Peschel, 
X. Wang, M. Kaulke, and K. Hallberg (Berlin: Springer, 
1999), pp. 91-125; M. A. Martm-Delgado, M. Roncaglia, 
and G. Sierra, Phys. Rev. B 64, 075117 (2001). 

8 N. Maeshima et ai, Phys. Rev. E 64 016705 (2001); Y. 
Nishio et ai, cond-mat/0401115, and references therein. 

9 F. Verstraete and J. I. Cirac, cond-mat/0407066; V. Murg, 
F. Verstaete, and J. I. Cirac, Phys. Rev. A 75, 033605 
(2007). 

10 A. Isacsson and O. F. Syljuasen, Phys. Rev. E 74, 026701 

(2006) ; V. Murg, F. Verstraete, and J. I. Cirac, Phys. Rev. 
B 79, 195119 (2009). 

11 Z. C. Gu, M. Levin, and X. G. Wen, Phys. Rev. B 78, 
205116 (2008); Z. C. Gu et ai, ibid 79, 085118 (2009). 

12 H. C. Jiang, Z. Y. Weng, and T. Xiang, Phys. Rev. Lett. 
101, 090603 (2008); Z. Y. Xie et al, arXiv:0809.0182v2. 

13 M. Levin and C. P. Nave, Phys. Rev. Lett. 99, 120601 

(2007) . 

14 M. Hinczewski and A. N. Berker, Phys. Rev. E 77, 011104 

(2008) ; M.-C. Chang and M.-F. Yang, Phys. Rev. B 79, 
104411 (2009). 



15 G. Vidal, Phys. Rev. Lett. 98, 070201 (2007); R. Orris and 
G. Vidal, Phys. Rev. B 78, 155117 (2008). 

16 M. Kohno and M. Takahashi, Phys. Rev. B 56, 3212 
(1997). 

17 S. Yunoki, Phys. Rev. B 65, 092402 (2002). 

18 It is known that the Hamiltonian in Eq. (|3]) is equivalent 
to a hard-core boson model under the mapping: a J = Sj + 

iff!-, a 3 = Sj-iS!j, and rij = aja,- = (SJ + 1/2), where a] is 
a creation operator of a hard-core boson at site j. Thus the 
same results are also found in terms of the boson language. 
For example, see G. G. Batrouni and R. T. Scalettar, Phys. 
Rev. Lett. 84, 1599 (2000); F. Hebert et al, Phys. Rev. B 
65, 014513 (2001). 

19 M. Suzuki, Phys. Lett. A 146, 319 (1990); J. Math Phys. 
(NY.) 32, 400 (1991). 

20 The Hamiltonain described in Eq. ((3| is equivalent to that 
of the conventional antiferromagnetic XXZ model by the 
following unitary transformation: rotating the spins along 
the 2 axis with an angle n (thus S™ —> —Sf and —* —S?) 
on the sublattice consisting of, say, even sites. 

21 S. Sachdev, Quantum Phase Transitions, (Cambridge Uni- 
versity Press, Cambridge, 1999). 

22 Within the N eel-ordered phase, one can show that 
S(Q)/N = |mf| 2 in the thermodynamic limit, where 
S(Q) = (1/N) J^i ■ S J )e iQ (ri - r -i ) is the spin structure 
factor. The QMC data for S(Q)/N in the thermodynamic 
limit can be read off horn Fig. 2(b) of Ref. [T3- 

23 For the spin-flopping states with the in-plane magnetiza- 
tion along the x direction, by using the relations between 
boson operators and the spin ones (see Ref. Il8h . we have 
po — \mZ\ 2 in the thermodynamic limit. Here po is the Bose 
condensate fraction in the hard-core boson model given by 
Po = (l/Y 2 )^ ij (a l t a,>. 



